1  Análisis de expresión génica diferencial (RNA-Seq)

El análisis de la expresión génica diferencial (RNA-Seq) es una de las técnicas bioinformáticas que más auge ha experimentado en la última década, ya que permite cuantificar la expresión de un gen a nivel celular y comparar estas expresiones entre distintos grupos experimentales para identificar asociaciones entre genes y factores experimentales (como enfermedades o tratamientos).

El proceso de un análisis genético diferencial tiene distintas etapas, que van desde la obtención del las células a estudiar, pasando por la secuenciación del RNA que contienen, hasta el análisis estadístico de estas frecuencias para determinar el nivel de expresión de cada gen y la comparación entre los grupos experimentales.

Etapas del análisis de expresión génica diferencial

En este tutorial nos centraremos en el análisis estadístico de las frecuencias de expresión de los genes, que normalmente requiere las siguientes etapas:

  1. Creación de la estructura de datos con las frecuencias de expresión génica.

  2. Preprocesamiento de datos.

  3. Análisis exploratorio de los datos.

  4. Análisis de expresión génica diferencial.

  5. Visualización de resultados.

Existen multitud de paquetes en el repositorio Bioconductor para realizar el análisis de expresión génica diferencial. En este tutorial veremos dos de los más usados: EdgeR y DesSeq2.

1.1 Análisis de expresión génica diferencial con edgeR

El paquete edgeR es uno de los principales paquetes usados para el análisis de expresión génica diferencial que incorpora funciones para todas las etapas del análisis estadístico. Está disponible en el repositorio Bioconductor.

A continuación se explica cómo realizar las distintas etapas del análisis estadístico con este paquete.

1.1.1 Estructura de datos

El análisis estadístico comienza siempre a partir de la matriz con las frecuencias o conteos de expresión de cada gen en el conjunto de muestras analizadas. Las filas de esta matriz corresponden a los genes observados y las columnas a las muestras analizadas.

Ejemplo 1.1 A continuación se muestran las primeras filas de una matriz de frecuencias de expresión génica. Cada casilla recoge el número de veces que el gen de la fila correspondiente se ha observado en la muestra de la columna correspondiente.

sample1 sample2 sample3 sample4 sample5 sample6
gene1 85 76 103 107 140 124
gene2 1 6 11 6 1 4
gene3 80 98 39 82 97 113
gene4 92 83 59 85 100 98
gene5 36 24 18 50 22 15
gene6 0 0 1 4 2 3

El paquete EdgeR utiliza la estructura de datos DGEList para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se utiliza la función

  • DGEList(frec, group = grupo): Crea una estructura del tipo DGEList con la matriz de frecuencias de expresión génica frec (con genes en filas y muestras en columnas) y el vector grupo con los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.
library(edgeR)
library(DEFormats)
frec <- simulateRnaSeqData()
grupo <- rep(c("A", "B"), each = 3)
dge <- DGEList(frec, group = grupo)
dge
An object of class "DGEList"
$counts
      sample1 sample2 sample3 sample4 sample5 sample6
gene1      85      76     103     107     140     124
gene2       1       6      11       6       1       4
gene3      80      98      39      82      97     113
gene4      92      83      59      85     100      98
gene5      36      24      18      50      22      15
995 more rows ...

$samples
        group lib.size norm.factors
sample1     A    42832            1
sample2     A    40511            1
sample3     A    39299            1
sample4     B    43451            1
sample5     B    40613            1
sample6     B    43662            1

Esta estructura de datos es una lista con dos atributos. El atributo counts contiene la matriz de frecuencias de expresión génica, y el atributo samples es un data frame con información sobre las muestras estudiadas. Cada fila de este data frame se corresponde con una columna de la matriz de frecuencias y contiene las siguientes columnas

Columna Descripción
group Grupo experimental al que pertenece la muestra.
lib.size tamaño de la librería (suma de frecuencias de la muestra).
norm.factors Factor de normalización.

La estructura de datos DGEList puede contener opcionalmente el atributo genes con anotaciones de los genes observados en las filas de la matriz de frecuencias.

Para convertir esta estructura de datos en una del tipo DESeqDataSet se puede utilizar la función as.DESeqDataSet del paquete DEFormats.

library(DESeq2)
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: 'matrixStats'
The following objects are masked from 'package:Biobase':

    anyMissing, rowMedians
The following object is masked from 'package:dplyr':

    count

Attaching package: 'MatrixGenerics'
The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars
The following object is masked from 'package:Biobase':

    rowMedians
dds <- as.DESeqDataSet(dge)
dds
class: DESeqDataSet 
dim: 1000 6 
metadata(1): version
assays(1): counts
rownames(1000): gene1 gene2 ... gene999 gene1000
rowData names(0):
colnames(6): sample1 sample2 ... sample5 sample6
colData names(3): group lib.size norm.factors

Ejemplo 1.2 Veamos cómo crear la estructura de datos correspondiente al estudio de Sheridan JM, et al. (2015) donde se obtuvieron datos de tres poblaciones de células (basal, progenitor luminal (LP) y luminal maduro (ML)) seleccionadas de las glándulas mamarias de ratones vírgenes hembra, cada una por triplicado. Los datos pueden obtenerse del repositorio Gene Expression Omnibus (GEO) mediante el identificador GSE63310.

En primer lugar cargamos la matriz de frecuencias.

frecuencias <- read.csv("datos/matriz-frecuencias-genes.csv", row.names = 1)
head(frecuencias)
          LP1 ML1   B1   B2  ML2 LP2  B3 ML3 LP3
497097      1   2  342  526    3   3 535   2   0
100503874   0   0    5    6    0   0   5   0   0
100038431   0   0    0    0    0   0   1   0   0
19888       0   1    0    0   17   2   0   1   0
20671       1   1   76   40   33  14  98  18   8
27395     431 771 1368 1268 1564 769 818 468 342

Después cargamos el diseño con los grupos experimentales.

diseño <- read.csv("datos/diseño-experimental.csv")
diseño
   Id Grupo Línea
1 LP1    LP  L004
2 ML1    ML  L004
3  B1 Basal  L004
4  B2 Basal  L006
5 ML2    ML  L006
6 LP2    LP  L006
7  B3 Basal  L006
8 ML3    ML  L008
9 LP3    LP  L008

A continuación creamos la estructura de datos DGEList.

# Creación de la estructura de datos DGEList
dge <- DGEList(counts = frecuencias, group = diseño$Grupo)
# Añadimos también la línea al diseño
dge$samples$Línea <- diseño$Línea
dge
An object of class "DGEList"
$counts
          LP1 ML1  B1  B2 ML2 LP2  B3 ML3 LP3
497097      1   2 342 526   3   3 535   2   0
100503874   0   0   5   6   0   0   5   0   0
100038431   0   0   0   0   0   0   1   0   0
19888       0   1   0   0  17   2   0   1   0
20671       1   1  76  40  33  14  98  18   8
27174 more rows ...

$samples
    group lib.size norm.factors Línea
LP1    LP 32863052            1  L004
ML1    ML 35335491            1  L004
B1  Basal 57160817            1  L004
B2  Basal 51368625            1  L006
ML2    ML 75795034            1  L006
LP2    LP 60517657            1  L006
B3  Basal 55086324            1  L006
ML3    ML 21311068            1  L008
LP3    LP 19958838            1  L008

A continuación anotamos los genes de las filas de la matriz de frecuencias. Para ello debe utilizarse un paquete específico para el genoma de la especie de donde proviene el RNA (Mus.muculus en este caso para el genoma del ratón).

library(Mus.musculus)
geneid <- rownames(dge)
genes <- select(Mus.musculus, keys = geneid, columns = c("SYMBOL", "TXCHROM"), keytype = "ENTREZID")
# Eliminamos duplicidad de algunos genes
genes <- genes[!duplicated(genes$ENTREZID),]
dge$genes <- genes
dge
An object of class "DGEList"
$counts
          LP1 ML1  B1  B2 ML2 LP2  B3 ML3 LP3
497097      1   2 342 526   3   3 535   2   0
100503874   0   0   5   6   0   0   5   0   0
100038431   0   0   0   0   0   0   1   0   0
19888       0   1   0   0  17   2   0   1   0
20671       1   1  76  40  33  14  98  18   8
27174 more rows ...

$samples
    group lib.size norm.factors Línea
LP1    LP 32863052            1  L004
ML1    ML 35335491            1  L004
B1  Basal 57160817            1  L004
B2  Basal 51368625            1  L006
ML2    ML 75795034            1  L006
LP2    LP 60517657            1  L006
B3  Basal 55086324            1  L006
ML3    ML 21311068            1  L008
LP3    LP 19958838            1  L008

$genes
   ENTREZID  SYMBOL TXCHROM
1    497097    Xkr4    chr1
2 100503874 Gm19938    <NA>
3 100038431 Gm10568    <NA>
4     19888     Rp1    chr1
5     20671   Sox17    chr1
27174 more rows ...

1.1.2 Preprocesamiento

Una vez preparada la estructura de datos, el siguiente paso es el preprocesamiento de datos. Normalmente comprende las siguientes tareas:

  1. Normalización de las frecuencias.
  2. Eliminación de los genes con poca expresión.
  3. Normalización de las distribuciones de frecuencias de expresión génicas.

1.1.2.1 Normalización de las frecuencias

El objetivo principal es normalizar las frecuencias de expresión génica para eliminar el efecto sobre las frecuencias de factores como la profundidad de secuenciado (a mayor profundidad de secuenciado mayores frecuencias) o el tamaño de los genes (a mayor tamaño mayores frecuencias). Generalmente se usan las siguientes transformaciones

  • Frecuencias por millón (CPM). Se divide cada frecuencia por el tamaño en millones de la librería de la muestra a la que pertenece. Se realiza con la función cpm(dge).

  • Logaritmo en base 2 de las frecuencias por millón (log2-CPM). Se calcula a partir de la anterior mediante la fórmula \(\log_2(CPM+\frac{2}{L})\), donde \(L\) es la longitud media de las librerías en millones. El término \(\frac{2}{L}\) que se añade permite calcular el logaritmo para frecuencias nulas (ya que el logaritmo de 0 no existe). Se realiza con la función cmp(dge, log = TRUE).

  • Lecturas por kilobase de transcripción por millón (RPKM): Se realiza con la función rpkm(dge, longitud_genes) pasando la longitud de los genes.

  • Fragmentos por kilobase de transcripción por millón (FPKM).

Las dos primeras no tienen en cuenta el tamaño de los genes, pero suelen usarse para la expresión génica diferencial donde el tamaño de los genes se supone constante en las distintas muestras.

Ejemplo 1.3 Siguiendo con el ejemplo anterior, calculamos las frecuencias por millón y el logaritmo de las frecuencias por millón.

cpm <- cpm(dge)
summary(cpm)
      LP1                 ML1                  B1                 B2          
 Min.   :    0.000   Min.   :    0.000   Min.   :   0.000   Min.   :   0.000  
 1st Qu.:    0.000   1st Qu.:    0.000   1st Qu.:   0.000   1st Qu.:   0.000  
 Median :    0.578   Median :    0.736   Median :   0.892   Median :   0.895  
 Mean   :   36.793   Mean   :   36.793   Mean   :  36.793   Mean   :  36.793  
 3rd Qu.:   19.536   3rd Qu.:   23.546   3rd Qu.:  24.221   3rd Qu.:  23.341  
 Max.   :27807.947   Max.   :11546.719   Max.   :7951.408   Max.   :7389.433  
      ML2                LP2                  B3                ML3          
 Min.   :   0.000   Min.   :    0.000   Min.   :   0.000   Min.   :   0.000  
 1st Qu.:   0.000   1st Qu.:    0.000   1st Qu.:   0.000   1st Qu.:   0.000  
 Median :   0.699   Median :    0.711   Median :   0.908   Median :   0.845  
 Mean   :  36.793   Mean   :   36.793   Mean   :  36.793   Mean   :  36.793  
 3rd Qu.:  23.827   3rd Qu.:   19.928   3rd Qu.:  21.439   3rd Qu.:  24.260  
 Max.   :7955.680   Max.   :29572.361   Max.   :9376.973   Max.   :7865.350  
      LP3           
 Min.   :    0.000  
 1st Qu.:    0.000  
 Median :    0.752  
 Mean   :   36.793  
 3rd Qu.:   21.594  
 Max.   :16500.710  
lcpm <- cpm(dge, log = TRUE)
summary(lcpm)
      LP1               ML1                B1                 B2         
 Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
 1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
 Median :-0.6847   Median :-0.3589   Median :-0.09513   Median :-0.0901  
 Mean   : 0.1714   Mean   : 0.3312   Mean   : 0.43559   Mean   : 0.4089  
 3rd Qu.: 4.2913   3rd Qu.: 4.5601   3rd Qu.: 4.60081   3rd Qu.: 4.5475  
 Max.   :14.7632   Max.   :13.4952   Max.   :12.95700   Max.   :12.8513  
      ML2               LP2                B3                ML3         
 Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
 1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
 Median :-0.4281   Median :-0.4064   Median :-0.07152   Median :-0.1704  
 Mean   : 0.3225   Mean   : 0.2529   Mean   : 0.40428   Mean   : 0.3708  
 3rd Qu.: 4.5772   3rd Qu.: 4.3199   3rd Qu.: 4.42513   3rd Qu.: 4.6031  
 Max.   :12.9578   Max.   :14.8520   Max.   :13.19491   Max.   :12.9413  
      LP3         
 Min.   :-4.5074  
 1st Qu.:-4.5074  
 Median :-0.3300  
 Mean   : 0.2749  
 3rd Qu.: 4.4355  
 Max.   :14.0102

1.1.2.2 Eliminación de genes con poca expresión

Aunque el objetivo del análisis de la expresión génica diferencial es detectar los genes que se expresan en un grupo experimental en comparación con los otros, cuando un gen no se expresa en ninguna de las muestras no tiene interés para el estudio y puede eliminarse.

Existen diferentes criterios para eliminar los genes con poca expresión. El paquete EdgeR incorpora la siguiente función

  • filterByExpr(dge): Devuelve un vector de booleanos con nombres correspondientes a los genes de la estructura DEGList dge. Por defecto devuelve TRUE para cualquier gen con una frecuencia mayor o igual que 10 en al menos un número de muestras igual al del menor grupo experimental.

Ejemplo 1.4 Veamos cuántos genes no tienen expresión en ninguna muestra del ejemplo anterior.

# Número de genes con expresión nula.
table(rowSums(dge$counts) == 0)  

FALSE  TRUE 
22026  5153

El siguiente gráfico muestra la distribución del logaritmo en base 2 de las frecuencias por millón.

# Definimos una función para dibujar la distribución del logaritmo de las frecuencias por millón.
dist_logcpm <- function(lcpm) {
lcpm |> 
    as.tibble()  |>
    pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
    ggplot(aes(x = LogCPM, color = Muestra)) +
    geom_density() +
    labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}

dist_logcpm(lcpm)
Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Como se observa ha una porcentaje significativo de genes que tiene una expresión muy baja (valores negativos).

A continuación eliminamos de la estructura de datos los genes con poca expresión.

filtro <- filterByExpr(dge)
dge <- dge[filtro, , keep.lib.sizes = F]

Y volvemos a dibujar la distribución del logaritmo en base 2 de las frecuencias por millón.

lcpm <- cpm(dge, log = TRUE)
dist_logcpm(lcpm)

Como se aprecia en el diagrama, el porcentaje de genes con baja expresión ha disminuido significativamente.

1.1.2.3 Normalización de las distribuciones de frecuencias de expresión génicas

Durante el proceso de secuenciación puede haber factores externos no biológicos que afecten a la expresión de los genes. Por ejemplo, la muestras del primer lote pueden tener mayores frecuencias que las del segundo lote. Como se supone que las muestras replicadas deben tener un rango y una distribución de frecuencias similares, en esta etapa se aplica otro procedimiento de normalización para garantizar que la distribución de frecuencias de cada muestra es similar. Para ello el paquete edgeR utiliza el método de las medias recortadas de los valores M por medio de la función

  • calcNormFactors(dge, method = "TMM"): Calcular los factores de escalado del tamaño de las librerías de cada muestra. Estos factores se guardan automáticamente en el data frame con la información de las muestras dge$samples.

Ejemplo 1.5 Siguiendo con el mismo ejemplo, calculamos los factores de escalado para cada muestra.

dge <- calcNormFactors(dge, method = "TMM")
dge$samples$norm.factors
[1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959
[8] 1.0861026 0.9839203

A continuación se muestran los diagramas de cajas de las distribuciones normalizadas tras aplicar los factores de escalado.

box_logcpm <- function(dge){
  diseño <- dge$samples
  diseño$Muestra <- row.names(diseño)
  cpm(dge, log = TRUE) |> 
    as.tibble() |>
    pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
    left_join(diseño, by = "Muestra") |> 
    ggplot(aes(x = Muestra, y = LogCPM, fill = group)) +
    geom_boxplot() +
    labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}

box_logcpm(dge)

Como se puede apreciar, después de la normalización, todas las muestras presentan una distribución de frecuencias similar.

1.1.3 Análisis exploratorio

Una vez preprocesados los datos comienza el análisis estadístico propiamente dicho. En una primera fase se suele realizar un análisis exploratorio de los datos que suele incluir los siguientes análisis:

  • Escalado multidimensional (análisis de componentes principales).

1.1.3.1 Escalado multidimensional

El escalado multidimensional mediante componente principales permite ver qué muestras son similares en cuando a la distribución de frecuencias de expresión génica. Los componentes principales son combinaciones lineales de los genes de la matriz de frecuencias que son ortogonales entre sí. El primer componente principal recoge la dimensión con mayor variabilidad de las frecuencias. El segundo recoge la dimensión don la mayor variabilidad no explicada por el primer componente principal y así sucesivamente. Normalmente los dos primeros componentes principales suelen recoger un porcentaje bastante alto de la variabilidad de las frecuencias. Al representar las muestras en estas dos dimensiones, las muestras más próximas entre sí, presentan una distribución de frecuencias de expresión génica similar. Para realizar esta representación se puede utilizar la siguiente función del paquete limma:

  • plotMDS(dge): Realiza un diagrama de componentes principales de la matriz de frecuencias de la estructura de datos dge.

Alternativamente, se pueden calcular los componentes principales mediante la función prcomp del paquete stats y luego usar la función autoplot del paquete ggfortify para dibujar el diagrama de componentes principales.

Ejemplo 1.6 A continuación se muestra cómo obtener el diagrama de componentes principales de la matriz de frecuencias de nuestro ejemplo.

plotMDS(dge, col = as.numeric(dge$samples$group), main = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")

library(ggfortify)
autoplot(prcomp(t(lcpm)), data = dge$samples, color = "group", loadings = F, loadings.label = F) +
labs(title = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")

Como se puede apreciar las muestras de cada grupo experimental aparecen agrupadas. La mayor diferencia (a lo largo del primer componente principal) se da entre el grupo basal y los otros dos grupos. Esto indica que cuando se realice el análisis de expresión génica diferencial habrá bastantes genes con diferencias de expresión significativa entre el grupo basal y los otros dos, mientras que no habrá tantos al comparar los grupos LP y ML.

Otra opción muy interesante es el paquete Glimma que permite dibujar un diagrama de componentes principales interactivo mediante la función

  • glMDSPlot(lcpm): Dibuja un diagrama de componentes principales interactivo de la matriz de frecuencias lcpm.

Ejemplo 1.7 A continuación se muestra cómo obtener el diagrama de componentes principales para nuestro ejemplo con el paquete Glimma.

library(Glimma)
#glMDSPlot(lcpm, groups = dge$samples[,c(2,5)])
dds <- DESeqDataSetFromMatrix(
  countData = dge$counts,
  colData = dge$samples,
  rowData = dge$genes,
  design = ~group
)
glimmaMDS(dds)

1.1.4 Análisis de expresión génica diferencial

La última etapa del análisis, y la más importante, es la detección de los genes con una expresión significativamente diferente entre los grupos experimentales. Para ello se suelen utilizar modelos lineales o modelos lineales generalizados. En general, la estimación de los parámetros del modelo ajustado depende de la distribución teórica usada para modelizar las frecuencias de expresión génica. El paquete limma, por ejemplo, realiza un ajuste de modelo lineal suponiendo que las distribuciones de las variables implicadas son normales, mientras que el paquete edgeR modeliza las frecuencias de expresión de los genes observadas mediante la distribución binomial negativa, que es mucho más realista.

En general en las distribuciones de frecuencias de expresión génica, se ha observado que la varianza no es independiente de la media. Los métodos que modelizan las frecuencias mediante el modelo binomial negativo asumen una relación cuadrática entre la varianza y media.

El siguiente paso es estimar la dispersión de las frecuencias de expresión génica para cada gen. Para ello edgeR utiliza el método de la máxima verosimilitud condicionada y ajustada por percentiles, mediante la función

  • estimateDisp(dge, diseño): Realiza una estimación de la dispersión de las frecuencias de expresión génica contenidas en la estructura del tipo DGEList dad en dge, de acuerdo al diseño experimental dado en diseño.

Ejemplo 1.8  

dge <- estimateDisp(dge)
Using classic mode.

Una vez ajustado el modelo Binomial Negativo y estimada la dispersión, solo queda aplicar el contraste de comparación de la expresión génica diferencial entre los grupos experimentales. Para ello edgeR utiliza el test exacto similar al test exacto de Fisher, pero adaptado a la distribución Binomial Negativa. Para ello se usa la función

  • exacTest(dge, pair=grupos): Realiza el contraste de comparación de la expresión génica a partir del las frecuencias de aparción de genes contenidas en la estructura del tipo DGEList dada en dge, entre los grupos experimentales dados en el vector grupos (solo puede contener dos grupos).

Ejemplo 1.9 En nuestro ejemplo realizamos tres contrastes para todos los posibles pares de grupos experimentales.

exact_B_LP <- exactTest(dge, pair= c("Basal", "LP"))
exact_B_ML <- exactTest(dge, pair= c("Basal", "ML"))
exact_LP_ML <- exactTest(dge, pair= c("LP", "ML"))

Tras realizar el contraste se puede utilizar la siguiente función para obtener un resumen con el número de genes subrexpresados y sobrexpresados significativamente.

  • summary(decideTests(test)): Devuelve una tabla con los genes subexpresados, sobreexpresados significativamente, así como lo que no presentan cambios significativos a partir de los resultados del contraste test.

Finalmente, para ver los genes con diferencias de expresión más estadísticamente significativas se puede utilizar la función

  • topTags(test): Devuelve un data frame con los genes con diferencias de expresión estadísticamente más significativas entre los grupos experimentales comparados en objeto de resultados del contraste test.

Ejemplo 1.10 Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo Basal y LP.

summary(decideTests(exact_B_LP))
       LP-Basal
Down       4864
NotSig     7203
Up         4557

Genes con diferencias más significativas entre el grupo Basal y el grupo LP.

topTags(exact_B_LP)
Comparison of groups:  LP-Basal 
       ENTREZID   SYMBOL TXCHROM     logFC   logCPM        PValue           FDR
67451     67451     Pkp2   chr16  5.661601 5.672839 2.184632e-164 3.631733e-160
19253     19253   Ptpn18    chr1  5.588760 5.344834 3.519102e-160 2.925078e-156
218518   218518 Marveld2   chr13  5.216899 6.054980 4.369286e-157 2.421167e-153
50722     50722    Dkkl1    chr7  6.296222 6.098528 6.043793e-156 2.511800e-152
227960   227960      Gca    chr2  5.516041 5.060151 6.024796e-144 2.003124e-140
242505   242505    Rasef    chr4  6.071991 6.704767 2.430700e-143 6.734659e-140
320007   320007    Sidt1   chr16  6.297688 4.897542 1.142473e-135 2.713210e-132
22249     22249   Unc13b    chr4  4.343846 6.600123 1.313443e-134 2.729335e-131
66871     66871    Cpne8   chr15 -4.685486 5.862087 4.580752e-129 8.461158e-126
269132   269132 Colgalt2    chr1 -4.989659 5.011534 4.878508e-126 8.110031e-123

Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo Basal y ML.

summary(decideTests(exact_B_ML))
       ML-Basal
Down       4773
NotSig     7011
Up         4840

Genes con diferencias más significativas entre el grupo Basal y el grupo ML.

topTags(exact_B_ML)
Comparison of groups:  ML-Basal 
       ENTREZID   SYMBOL TXCHROM     logFC   logCPM        PValue           FDR
50722     50722    Dkkl1    chr7  6.728394 6.098528 4.166088e-173 6.925705e-169
242505   242505    Rasef    chr4  6.678297 6.704767 2.371625e-165 1.971295e-161
21804     21804  Tgfb1i1    chr7 -6.120166 6.012779 9.151423e-156 5.071108e-152
22249     22249   Unc13b    chr4  4.561206 6.600123 3.599522e-150 1.495961e-146
20661     20661    Sort1    chr3  4.948551 7.722875 8.615672e-144 2.864539e-140
12521     12521     Cd82    chr2  4.666737 8.007651 8.771904e-143 2.430402e-139
67451     67451     Pkp2   chr16  5.100406 5.672839 7.763029e-140 1.843608e-136
218518   218518 Marveld2   chr13  4.815118 6.054980 7.247449e-139 1.506020e-135
78896     78896    Ecrg4    chr1 -6.587115 5.670792 9.526597e-139 1.759668e-135
214968   214968   Sema6d    chr2 -5.935614 6.101960 1.913875e-125 3.181626e-122

Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo LP y ML.

summary(decideTests(exact_LP_ML))
       ML-LP
Down    2432
NotSig 11226
Up      2966

Genes con diferencias más significativas entre el grupo LP y el grupo ML.

topTags(exact_LP_ML)
Comparison of groups:  ML-LP 
       ENTREZID  SYMBOL TXCHROM     logFC   logCPM       PValue          FDR
94212     94212    Pag1    chr3 -2.616277 5.916572 9.314994e-31 1.548525e-26
214968   214968  Sema6d    chr2 -2.378756 6.101960 3.482661e-27 2.894788e-23
114479   114479  Slc5a5    chr8  2.879278 7.770145 1.246288e-26 6.906100e-23
207592   207592 Tbc1d16   chr11  2.308784 5.659683 6.219699e-26 2.584907e-22
13617     13617   Ednra    chr8 -3.187364 4.014091 1.557011e-25 5.176750e-22
231605   231605  Galnt9    chr5  3.289480 2.482287 2.766785e-24 7.665838e-21
12614     12614  Celsr1   chr15  2.826770 5.693953 1.525133e-23 3.621973e-20
14778     14778    Gpx3   chr11  3.045448 9.996423 4.769484e-23 9.910988e-20
107895   107895   Mgat5    chr1  2.025411 5.535760 9.509983e-23 1.756600e-19
320127   320127    Dgki    chr6  2.415427 4.951953 1.261948e-22 2.097863e-19

Otra alternativa, para diseños experimentales que incluye más de un factor, es usar modelos lineales generalizados (GLM). El primer paso es definir la matriz del diseño del experimento, que incluye las variables que definen grupos experimentales. Para ello se utiliza la función

  • model.matrix(formula-diseño): Construye una matriz con el diseño experimental de acuerdo a la fórmula dada en formula-diseño. Esta fórmula es similar a la fórmula que que utiliza en un ANOVA para definir el diseño experimental, donde las variables se añaden con el operador +, y se hacen interactuar con el operador *.

Ejemplo 1.11 Continuando con el mismo ejemplo, definimos la matriz de diseño del modelo.

diseño <- model.matrix(~ 0 + dge$samples$group + dge$samples$Línea)
colnames(diseño) <- gsub("dge\\$samples\\$group", "", colnames(diseño))
colnames(diseño) <- gsub("dge\\$samples\\$", "", colnames(diseño))
diseño
  Basal LP ML LíneaL006 LíneaL008
1     0  1  0         0         0
2     0  0  1         0         0
3     1  0  0         0         0
4     1  0  0         1         0
5     0  0  1         1         0
6     0  1  0         1         0
7     1  0  0         1         0
8     0  0  1         0         1
9     0  1  0         0         1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"

attr(,"contrasts")$`dge$samples$Línea`
[1] "contr.treatment"

Al igual que antes, antes de realizar el contraste hay que estimar la dispersión conjunta y para cada gen, pero ahora hay que introducir el diseño como un parámetro de la función estimateDisp.

Ejemplo 1.12 Realizamos la estimación de la dispersión para nuestro ejemplo.

dge <- estimateDisp(dge, diseño)

Una vez estimada la dispersión y ajustados los modelos generalizados binomiales negativos, ya se puede realizar el contraste de comparación de expresión génica. Para ello se utilizan las siguientes funciones

  • glmQLFit(dge, diseño): Realiza el ajuste del modelo generalizado de comparación de expresión génica mediante el test F de cuasi verosimilitud, con las frecuencias de expresión de la estructura del tipo DGEList dada en dge y de acuerdo al diseño experimental indicado en diseño.

  • makeContrast(grupo1 - grupo2, levels=diseño): Define los grupos grupo1 y grupo2 a comparar en el contraste pares para la diferencia en la expresión génica.

  • glmQLFTest(modelo-ajustado, contrast=contraste). Realiza el contraste de comparación por pares a partir del modelo ajustado modelo-ajustado. El modelo ajustado tiene varios coeficientes dependiendo del número grupos experimentales, el primero corresponde al ajuste base para el grupo control, y los sucesivos a las diferencias de los otros grupos experimentales con el control.

  • glmTreat(modelo-ajustado, contrast=contraste, lfc = n). Realiza el contraste de comparación por pares similar al de la función anterior, pero descartando los genes con un logaritmo en base dos de la razón de cambio (logFC) menor que el valor dado en el parámetro lfc.

Ejemplo 1.13 Realizamos el ajuste del modelo linear generalizado para nuestro ejemplo.

modelo <- glmQLFit(dge, diseño)

Y, a continuación, los contrastes por pares. En primer lugar, comparamos Basal con LP y mostramos los genes más diferenciados.

contraste <- makeContrasts(Basal - LP, levels = diseño) 
glmQL_B_LP <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_B_LP))
       1*Basal -1*LP
Down            4523
NotSig          7192
Up              4909
topTags(glmQL_B_LP)
Coefficient:  1*Basal -1*LP 
       ENTREZID   SYMBOL TXCHROM     logFC   logCPM        F       PValue
242505   242505    Rasef    chr4 -5.941271 6.704195 1863.844 5.629045e-11
67451     67451     Pkp2   chr16 -5.745868 5.671932 1845.036 5.867451e-11
19253     19253   Ptpn18    chr1 -5.655466 5.343492 1616.133 1.008466e-10
53624     53624    Cldn7   chr11 -5.527387 7.529234 1327.691 2.251566e-10
14275     14275    Folr1    chr7 -6.925474 5.510509 1313.864 2.349910e-10
218518   218518 Marveld2   chr13 -5.153654 6.054055 1300.770 2.448002e-10
70350     70350    Basp1   chr15 -6.086771 6.624263 1223.255 3.145873e-10
228543   228543     Rhov    chr2 -6.263337 6.939962 1112.552 4.632731e-10
70737     70737      Cgn    chr3 -5.466154 6.644032 1075.182 5.325457e-10
320007   320007    Sidt1   chr16 -6.345681 4.895659 1071.587 5.398676e-10
                FDR
242505 4.877025e-07
67451  4.877025e-07
19253  5.588246e-07
53624  6.258841e-07
14275  6.258841e-07
218518 6.258841e-07
70350  6.258841e-07
228543 6.258841e-07
70737  6.258841e-07
320007 6.258841e-07

A continuación comparamos Basal con ML.

contraste <- makeContrasts(Basal - ML, levels = diseño) 
glmQL_B_ML <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_B_LP))
       1*Basal -1*LP
Down            4523
NotSig          7192
Up              4909
topTags(glmQL_B_ML)
Coefficient:  1*Basal -1*ML 
       ENTREZID   SYMBOL TXCHROM     logFC   logCPM        F       PValue
242505   242505    Rasef    chr4 -6.539106 6.704195 2151.638 3.128356e-11
71740     71740  Nectin4    chr1 -5.586420 6.447456 1566.715 1.144948e-10
67451     67451     Pkp2   chr16 -5.178572 5.671932 1566.704 1.144981e-10
53624     53624    Cldn7   chr11 -5.504748 7.529234 1333.741 2.210153e-10
78896     78896    Ecrg4    chr1  6.456502 5.672665 1256.914 2.815928e-10
19253     19253   Ptpn18    chr1 -4.594745 5.343492 1155.825 3.964889e-10
50722     50722    Dkkl1    chr7 -6.778277 6.097962 1152.189 4.016182e-10
12521     12521     Cd82    chr2 -4.687198 8.007496 1151.233 4.029809e-10
218518   218518 Marveld2   chr13 -4.755942 6.054055 1146.641 4.096058e-10
108153   108153  Adamts7    chr9  4.697586 5.748726 1135.707 4.259335e-10
                FDR
242505 5.200579e-07
71740  6.188524e-07
67451  6.188524e-07
53624  6.188524e-07
78896  6.188524e-07
19253  6.188524e-07
50722  6.188524e-07
12521  6.188524e-07
218518 6.188524e-07
108153 6.188524e-07

Y ahora comparamos LP con ML.

contraste <- makeContrasts(LP - ML, levels = diseño) 
glmQL_LP_ML <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_B_LP))
       1*Basal -1*LP
Down            4523
NotSig          7192
Up              4909
topTags(glmQL_LP_ML)
Coefficient:  1*LP -1*ML 
      ENTREZID  SYMBOL TXCHROM    logFC   logCPM        F       PValue
11815    11815    Apod   chr16 4.283581 6.321010 486.9859 1.330736e-08
20319    20319   Sfrp2    chr3 2.753550 4.598845 314.7591 7.723980e-08
15212    15212    Hexb   chr13 2.909540 5.986395 286.4532 1.126803e-07
13132    13132    Dab2   chr15 2.557433 5.184149 271.5224 1.395871e-07
14962    14962     Cfb   chr17 1.803855 4.762677 253.8898 1.825182e-07
18552    18552   Pcsk5   chr19 2.158846 4.236652 238.4773 2.342784e-07
12424    12424     Cck    chr9 4.313973 4.761400 236.2864 2.430493e-07
74365    74365  Lonrf3    chrX 2.245457 4.517287 226.4194 2.880104e-07
73341    73341 Arhgef6    chrX 2.346351 7.760006 219.2708 3.271778e-07
18858    18858   Pmp22   chr11 1.724631 5.964828 204.3723 4.325480e-07
               FDR
11815 0.0002212215
20319 0.0005622194
15212 0.0005622194
13132 0.0005622194
14962 0.0005622194
18552 0.0005622194
12424 0.0005622194
74365 0.0005622194
73341 0.0005622194
18858 0.0005622194

Por último buscamos los genes con una expresión diferente en los tres grupos experimentales.

glmQL_all <- glmQLFTest(modelo, coef = 2:3)
topTags(glmQL_all)
Coefficient:  LP ML 
       ENTREZID  SYMBOL TXCHROM  logFC.LP  logFC.ML   logCPM        F
68598     68598  Dnajc8    chr4 -14.75481 -14.52690 5.503444 3917.603
78658     78658  Ncapd3    chr9 -15.14518 -15.10475 5.096621 3893.545
76251     76251 Ercc6l2   chr13 -14.63388 -14.74224 5.220876 3886.831
67444     67444   Ilkap    chr1 -14.43347 -14.33013 5.399407 3841.724
58859     58859  Efemp2   chr19 -14.80616 -14.65589 5.410953 3840.347
24045     24045  Scamp3    chr3 -15.06093 -14.58439 5.202811 3765.952
21371     21371    Tbca   chr13 -14.99664 -15.05218 5.185991 3765.800
57443     57443   Fbxo3    chr2 -14.01332 -14.04999 5.861901 3765.551
80517     80517 Herpud2    chr9 -14.26672 -14.05908 5.853877 3756.947
231464   231464  Cnot6l    chr5 -14.28498 -14.01301 5.835142 3750.520
             PValue          FDR
68598  5.798476e-13 3.280446e-11
78658  5.946855e-13 3.280446e-11
76251  5.989106e-13 3.280446e-11
67444  6.282843e-13 3.280446e-11
58859  6.292087e-13 3.280446e-11
24045  6.817772e-13 3.280446e-11
21371  6.818903e-13 3.280446e-11
57443  6.820748e-13 3.280446e-11
80517  6.885051e-13 3.280446e-11
231464 6.933572e-13 3.280446e-11

1.1.5 Visualización de resultados

Existen diferentes diagramas para representar los resultados del análisis génico diferencial. En esta sección presentamos dos de los diagramas más utlizados, el diagrama MD y el diagrama de volcán.

1.1.5.1 Diagrama MD

Este diagrama consiste en un mapa de puntos donde cada punto corresponde a un gen. El eje X representa la media del logaritmo de las frecuencias por millón (lcpm) y el eje Y representa el logaritmo en base 2 de la razón de cambio (logFC). Habitualmente, los genes subexpresados significativamente se dibujan en color azul, mientras que los sobreexpresados se dibujan en color rojo.

Para dibujar un diagrama MD con edgeR se utiliza la función

  • plotMD(contraste): Dibuja un diagrama MD a partir de los resultados del contraste de comparación de expresión génica dado en contraste.

Ejemplo 1.14 Dibujamos primero el diagrama MD para la comparación de los grupos Basal y LP.

plotMD(glmQL_B_LP)

A continuación para la comparación de los grupos Basal y ML.

plotMD(glmQL_B_ML)

Y finalmente para la comparación de los grupos LP y ML.

plotMD(glmQL_LP_ML)

1.1.6 Diagrama de volcán

Un diagrama de volcán es una representación gráfica utilizada en Bioinformática para visualizar los resultados del análisis de expresión génica diferencial u otros tipos de análisis de datos ómicos de alto rendimiento, como Proteómica o Metabolómica.

En un diagrama de volcán, cada punto de datos representa un gen (o una proteína/metabolito) del conjunto de datos. El eje x muestra el cambio logarítmico o tamaño del efecto, que mide la magnitud del cambio en la expresión entre dos condiciones (por ejemplo, tratamiento vs. control). El eje y muestra la significación estadística, a menudo representada como el logaritmo negativo del valor p. Un valor p más pequeño indica una mayor significación estadística.

A menudo, los puntos en el gráfico de volcán están coloreados o resaltados según su significación estadística y cambio en la expresión. Por lo general, los genes significativamente sobreexpresados se representan en rojo, los genes significativamente subexpresados en azul y los genes no significativamente expresados en gris o negro. Cuanto más significativos sean estadísticamente y mayor sea el cambio en la expresión, más alejados estarán los puntos del centro del gráfico.

edgR no incorpora una función para realizar este tipo de diagramas, así que utilizaremos la siguiente función del paquete Glimma

  • glimmaVolcano(modelo, dge = dge): Realiza diagrama de volcán del modelo ajustado fit sobre la estructura de datos del tipo DGEList dada en dge.

En el capítulo Diagramas de Volcán se explica con más detalle su creación con el paquete ggplot2.

Ejemplo 1.15 Dibujamos primero el diagrama de volcán para la comparación de los grupos Basal y LP.

library(Glimma)
glimmaVolcano(glmQL_B_LP, dge = dge)

A continuación para la comparación de los grupos Basal y ML.

library(Glimma)
glimmaVolcano(glmQL_B_ML, dge = dge)

Y finalmente para la comparación de los grupos LP y ML.

library(Glimma)
glimmaVolcano(glmQL_LP_ML, dge = dge)

1.2 Análisis de expresión génica diferencial con DESeq2

El paquete DESeq2 es otro de los paquetes más utilizados para el análisis de la expresión génica diferencial, que al igual que edgeR incorpora procedimientos para todas las etapas del análisis. También está disponible en el repositorio Bioconductor.

1.2.1 Estructura de Datos

El paquete DESeq2 utiliza la estructura de datos DESeqDataSet para guardar la matriz de frecuencias de expresión génica. Esta estructura se deriva, a su vez, de la estructura SingleCellExperiment, que además de la matriz de frecuencias de expresión de los genes incluye también los grupos experimentales a los que pertenecen las muestras estudiadas y otra información que se va generando a medida que progresa el análisis.

El objeto SingleCellExperiment

El paquete SingleCellExperiment de Bioconductor define esta estructura de datos, y puede instalarse mediante el siguiente código

BiocManager::install('SingleCellExperiment')

Para crear la estructura de datos DESeqDataSet se utiliza la siguiente función

  • DESeqDataSetFromMatrix(countData = frec, colData = grupo, design = diseño): Crea una estructura del tipo DESeqDataSet con la matriz de frecuencias de expresión génica frec (con genes en filas y muestras en columnas), y el data frame grupo con los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias y la columna diseño del data frame grupo que contiene los grupos experimentales a comparar.
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
  countData = dge$counts,
  colData = dge$samples,
  rowData = dge$genes,
  design = ~group
)
dds
class: DESeqDataSet 
dim: 16624 9 
metadata(1): version
assays(1): counts
rownames(16624): 497097 20671 ... 100861924 170942
rowData names(3): ENTREZID SYMBOL TXCHROM
colnames(9): LP1 ML1 ... ML3 LP3
colData names(4): group lib.size norm.factors Línea

1.2.2 Preprocesamiento

El paquete DESeq2 realiza el preprocesamiento de maneare automática cuando se lanza el procedimiento para el análisis de expresión génica diferencial.

1.2.3 Análisis exploratorio

1.2.4 Escalado multidimensional

Al igual que antes, el principal análisis exploratorio es el escalado multidimensional. Para realizarlo podemos usar de nuevo la siguiente función del plaqute Glimma.

  • glimmaMDS(dds): Realiza un diagrama interactivo de componenetes principales a partir de la estructura de datos dds.

Ejemplo 1.16 A continuación se muestra cómo obtener el diagrama de componentes principales para nuestro ejemplo con el paquete Glimma.

library(Glimma)
glimmaMDS(dds)

1.2.5 Análisis de expresión génica diferencial

Para realizar el análisis de expresión génica diferencial el paquete DESeq2 utiliza la siguiente función

  • DESeq(dds): Realiza el preprocesamiento de datos y el ajuste del modelo para el análisis de expresión génica diferencial de la estructura de datos dds.

Ejemplo 1.17 A continuación se muestra cómo realizar el análisis de expresión génica diferencial para nuestro ejemplo.

dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
dds
class: DESeqDataSet 
dim: 16624 9 
metadata(1): version
assays(4): counts mu H cooks
rownames(16624): 497097 20671 ... 100861924 170942
rowData names(29): ENTREZID SYMBOL ... deviance maxCooks
colnames(9): LP1 ML1 ... ML3 LP3
colData names(5): group lib.size norm.factors Línea sizeFactor

1.2.6 Visualización de datos

1.2.6.1 Diagrama MA

Ejemplo 1.18  

glimmaMA(dds)
15 of 16624 genes were filtered out in DESeq2 tests

1.2.6.2 Diagrama de volcán

Ejemplo 1.19  

glimmaVolcano(dds)

1.3 Referencias